Event-by-event anisotropic flow in heavy-ion collisions 
from combined Yang-Mills and viscous fluid dynamics 
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Anisotropic flow coefficients V1-V5 in heavy ion collisions are computed by combining a classi- 
cal Yang-Mills description of the early time glasma flow with the subsequent relativistic viscous 
hydrodynamic evolution of matter through the quark-gluon plasma and hadron gas phases. The 
glasma dynamics, as realized in the IP-Glasma model, takes into account event-by-event geometric 
fluctuations in nucleon positions and intrinsic sub-nucleon scale color charge fluctuations; the pre- 
equilibrium flow of matter is then matched to the MUSIC algorithm describing viscous hydrodynamic 
flow and particle production at freeze-out. The IP-Glasma+MUSlC model describes well both trans- 
verse momentum dependent and integrated v n data measured at the Large Hadron Collider (LHC) 
and the Relativistic Heavy Ion Collider (RHIC). The model also reproduces the event-by-event dis- 
tributions of i>2, ^3 and V4 measured by the ATLAS collaboration. The implications of our results 
for better understanding of the dynamics of the glasma as well as for the extraction of transport 
properties of the quark-gluon plasma are outlined. 



Heavy ion collisions at the Relativistic Heavy Ion Col- 
lider (RHIC) and the Large Hadron Collider (LHC) 
uniquely allow for systematic exploration of the high tem- 
perature many-body dynamics of a non-Abclian quan- 
tum field theory. Particularly intriguing is the prospect 
of disentangling the non-equilibrium strongly correlated 
dynamics of the early time glasma regime from those of 
late stage nearly equilibrated quark-gluon plasma and 
hadron gas phases by measurements of anisotropic flow 
harmonics v n at both RHIC [1, 2] and LHC [3-5]. 

An excellent candidate for providing initial conditions 
for systematic flow studies is the IP-Glasma model de- 
scribed in detail in Refs. [6, 7]. It combines the IP-Sat 
(Impact Parameter Saturation Model) model [8, 9] of 
high energy nucleon (and nuclear) wavefunctions with the 
classical Yang-Mills (CYM) dynamics of the glasma fields 
produced in a heavy-ion collision [10-13]. We note that 
the IP-Sat model provides a good description of small x 
HERA deeply inelastic scattering (DIS) data off protons 
and fixed target nuclear DIS data [14]. Prior implemen- 
tation of the IP-Sat model in proton-proton and nucleus- 
nucleus collisions at the LHC using a fcj_-factorized ex- 
pression approximating CYM dynamics was shown to 
give good agreement with bulk features of data [15]. The 
upcoming p+Pb run at the LHC should provide further 
constraints on the dynamics of the IP-Glasma model, in 
particular the energy dependence of Q s . 

In this letter, we couple the IP-Glasma model of the 
classical early time evolution of boost-invariant config- 
urations of gluon fields to a relativistic hydrodynamic 
description of the system, using the energy density and 
flow velocity in the transverse plane at the switching time 
''"switch ~ l/<3s as input. The hydrodynamic evolution in 
each event is described by MUSIC [16-19], a 3+1 dimen- 
sional relativistic viscous hydrodynamic simulation [20] 
that uses the Kurganov-Tadmor algorithm [21]. While 



this matching of glasma dynamics to viscous hydrody- 
namics is a significant improvement relative to previ- 
ously employed initial conditions for heavy ion collisions, 
early stage dynamics is not fully included. Most notably, 
the hydrodynamic viscous tensor LP" is too large to be 
described sclf-consistcntly by a gradient expansion. In- 
stabilities triggered by quantum fluctuations, and subse- 
quent strong scattering of over-occupied fields, may lead 
to rapid quenching of YL^ to reasonable values justify- 
ing the use of viscous hydrodynamics already at early 
times. In this letter, we will assume such an efficient 
mechanism to be at work and set the initial value of 
to zero. Recent progress in computing early-time 
quantum fluctuations will help eliminate this systematic 
uncertainty [22-26]. 

When we switch from the CYM description to hydro- 
dynamics we construct the fluid's initial energy momen- 
tum tensor T^ id = (e + V)u^u v - Vcf u + 11^ from the 
energy density in the fluid's rest frame e, the flow veloc- 
ity w^, and, using an equation of state, the local pressure 
V at each transverse position, e and u M are obtained 
by solving u^T^ym ~ euh ' 1 using the fact that m' 1 is a 
time- like eigenvector of Tq YM and satisfies u 2 = 1. 

Other important details of our analysis are as follows. 
Unless otherwise noted, T switc h = 0.2fm/c. We employ 
the s95p-PCE equation of state, obtained from fits to lat- 
tice QCD results and a hadron resonance gas model [27], 
with partial chemical equilibrium (PCE) setting in be- 
low a temperature Tpce = 150 MeV. Kinetic freeze-out 
occurs at T-po — 120 McV. At this temperature, we im- 
plement the Cooper-Frye prescription [28] for computing 
particle spectra. Unless otherwise noted, shown results 
include decays from resonances of masses up to 1.3 GeV. 

A novel feature of our study is the determination of 
centrality classes using the multiplicity distribution of 
gluons much alike the procedure followed by the heavy 
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FIG. 1. (Color online) Gluon multiplicity distribution in the 
IP-Glasma model. 



10~ 



ALICE 

0-5% 



IP-Glasma+MUSIC - 

7i + + n m 

K + + K" v 

P + P ° 




FIG. 2. (Color online) Identified particle transverse momen- 
tum spectra including all resonances up to 2 GeV compared 
to experimental data from the ALICE collaboration [31]. 



ion experiments [29] . The gluon multiplicity distribution 
is shown in Fig. 1. Centrality classes arc determined from 
the fraction of the integral over this distribution, begin- 
ning with integrating from the right. As a consequence 
of implementing this centrality selection, we properly ac- 
count for impact parameter and multiplicity fluctuations. 

Because entropy is produced during the viscous hydro- 
dynamic evolution, we need to adjust the normalization 
of the initial energy density commensurately to describe 
the final particle spectra [30]. The obtained py-spectra 
of pions, kaons, and protons are shown for 0-5% central 
collisions at = 2.76 TeV/nucleon, using rj/s = 0.2, 
in Fig. 2, and compared to data from ALICE [31]. The 
results are for averages over only 20 events in this case, 
but statistical errors are smaller than the line width for 
the spectra. Overall, the agreement with experimental 
data is good. However, soft pions at pr < 300 MeV are 
underestimated . 

We determine v\ to v§ in every event by first deter- 
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FIG. 3. (Color online) Root-mean-square anisotropic flow co- 
efficients (v 2 ,) 1 ^ 2 as a function of transverse momentum, com- 
pared to experimental data by the ATLAS collaboration using 
the event plane (EP) method [4] (points). 200 events. Bands 
indicate statistical errors. Experimental error bars are smaller 
than the size of the points. 

mining the exact event plane [32] 



1 (sm(n(j))) 

ip n = - arctan — — 

n {cos(n<p)) 



(1) 



and then computing 



Vn(pr) = (cos(n(</> - tp n ))) 

_ J dcj)f(p ± ,(j)) cos(n((j) ~ ip n )) 

r#/(p ± ,0) 



(2) 



where f(p±, <p) are the thermal distribution functions ob- 
tained in the Cooper-Frye approach (with additional con- 
tributions from resonance decays). 

We first present the root-mean-square (rms) v n (pT) for 
10 — 20% central collisions and compare to experimental 
data from the ATLAS collaboration [4] in Fig. 3. Agree- 
ment for V2-V5 is excellent. We note that the v n from 
the experimental event plane method do not exactly cor- 
respond to the rms values, but lie somewhere between 
the mean and the rms values. In this regard, a better 
comparison is the p^-integrated rms v n to the ALICE 
t> n {2} results-which correspond to the rms values. Ex- 
cellent agreement over the whole studied centrality range 
is achieved for the experimentally available 1)2, v% and v±, 
as shown in Fig. 4. 

We studied the effect of initial transverse flow included 
in our framework by also computing v n (pr) with u M set 
to zero at time r sw i tc h- The effect on hadron anisotropic 
flow turns out to be extremely weak - results agree within 
statistical errors. Because photons are produced early 
on in the collision, we expect a greater effect on photon 
anisotropic flow; this will be examined in a subsequent 
work. We emphasize that prc-cquilibrium dynamics that 
is not fully accounted for may still influence the amount 
of initial transverse flow. 
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FIG. 4. (Color online) Root-mean-square anisotropic flow co- 
efficients {^) 1//2 , computed as a function of centrality, com- 
pared to experimental data of v n {2}, n £ {2,3,4}, by the 
ALICE collaboration [3] (points). Results are for 200 events 
per centrality with bands indicating statistical errors. 



0.5 1 1.5 2 

Pt [GeV] 

FIG. 6. (Color online) Comparison of v n (pT) using con- 
stant rj/s = 0.2 and a temperature dependent rj/s(T) as 
parametrized in [33]. Experimental data by the ATLAS col- 
laboration using the event-plane (EP) method [4] (points). 
Bands indicate statistical errors. 
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FIG. 5. (Color online) Comparison of u„(pt) using two dif- 
ferent switching times r 8W it c h = 0.2 fm/c (wide), and 0.4 fm/c 
(narrow). Experimental data by the ATLAS collaboration us- 
ing the event-plane (EP) method [4] (points). Bands indicate 
statistical errors. 

The effect of changing the switching time from 
^switch = 0.2 fm/c to T switch = 0.4 fm/c is shown in Fig. 5. 
Results agree within statistical errors, but tend to be 
slightly lower for the later switching time. The nonlinear 
interactions of classical fields become weaker as the sys- 
tem expands and therefore Yang-Mills dynamics is less 
effective than hydrodynamics in building up flow at late 
times. Yet it is reassuring that there is a window in time 
where both descriptions produce equivalent results. 

Because a constant rj/s is at best a rough effective 
measure of the evolving shear viscosity to entropy den- 
sity ratio, we present results for a parametrized temper- 
ature dependent rj/s, following [33]. We use the same 
paramctrization (HH-HQ) as in [33, 34] with a minimum 
of rj/s{T) = 0.08 at T = T tI = 180 MeV. The result, 



compared to rj/s = 0.2 is shown for 20 — 30% central col- 
lisions in Fig. 6. The results are indistinguishable when 
studying just one collision energy. The insensitivity of 
our results to two very different functional forms may 
suggest that a very large fraction of the magnitude of 
the flow coefficients is built up at later times when rj/s 
is very small. Also, since second order viscous hydrody- 
namics breaks down when IT A " y is comparable to the ideal 
terms, our framework may be inadequate for large values 
of rj/s. 

At top RH1C energy, as shown in Fig. 7, the experi- 
mental data from STAR [35] and PHENIX [1] is well de- 
scribed when using a constant rj/s = 0.12, which is about 
40 % smaller than the value at LHC. A larger effective rj/s 
at LHC than at RHIC was also found in [36]. The tem- 
perature dependent rj/s(T) used to describe LHC data 
works well for low-px RHIC data, but underestimates 
V2(pr) and vs(j)t) forpT > 1 GeV. The paramctrizations 
of rj/s(T) in the literature are not definitive and signif- 
icant improvements are necessary. Our studies suggest 
great potential for extracting the temperature dependent 
properties of QCD transport coefficients by performing 
complementary experiments extracting flow harmonics at 
both RHIC and LHC. 

In Fig. 8 we present results for v\(pt) compared to ex- 
perimental data from ALICE [37] , extracted in [39] , and 
from ATLAS [38]. v\(pr) cannot be positive definite be- 
cause momentum conservation requires (vi(pt)pt) = 0. 
There is a disagreement between the experimental results 
(discussed in [38] ) and between theory and experiment at 
LHC. On the other hand, v\(pt) at RHIC is very well re- 
produced (see Fig. 7). One possible explanation for the 
data crossing v\(pt) = at a lower pt than the calcu- 
lation at LHC could be the underestimation of the pion 
PT-spectra at very low px - see Fig. 2. However, this is 
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FIG. 7. (Color online) Comparison of v„(pt) at RHIC using 
constant r]/s = 0.12 and a temperature dependent rj/s(T) as 
parametrized in [33]. Experimental data by the PHENIX [1] 
(open symbols) and STAR [35] (preliminary, filled symbols) 
collaborations. Bands indicate statistical errors. 
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FIG. 8. (Color online) Vi(pr) compared to experimental data 
from the ALICE [37] and ATLAS [38] collaborations. 

not necessarily the only explanation. In fact, for RHIC 
energies, calculated pion spectra also underestimate the 
data for pt < 300 MeV but vi(px) is well reproduced. 

We present event-by-event distributions of i>2, «3, and 
V4 compared to results from the ATLAS collaboration 
[40, 41] in Fig. 9. We chose 20-25% central events be- 
cause eccentricity distributions from neither MC-Glauber 
nor MC-KLN models agree with the experimental data 
in this bin [41]. To compare data with the distribution 
of initial eccentricities [42] from the IP-Glasma model 
and the final v n distributions after hydrodynamic evolu- 
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FIG. 9. (Color online) Scaled distributions of V2, «3, and v± 
(from top to bottom) compared to experimental data from 
the ATLAS collaboration [40, 41]. 1300 events. Bands are 
systematic experimental errors. 



tion, we scaled the distributions by their respective mean 
value. We find that the initial eccentricity distributions 
are a good approximation to the distribution of experi- 
mental v n . Only for v 4 (and less so for V2) the large v n 
end of the experimental distribution is much better de- 
scribed by the hydrodynamic v n distribution than the e n 
distribution. This can be explained by non-linear mode 
coupling becoming important for large values of i>2 and 
v 4 . 

In summary, we have shown that the IP- 
Glasma+MUSiC model gives very good agreement 
to multiplicity and flow distributions at RHIC and LHC. 
By including properly sub-nucleon scale color charge 
fluctuations and their resulting early time CYM dynam- 
ics, this model significantly extends previous studies in 
the literature [19, 36, 43-47]. Omitted in all studies 
including ours is the stated dynamics of instabilities and 
strong scattering in over-occupied classical fields that 
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can drive the system to isotropy and generate substan- 
tial flow well prior to thermalization. Ongoing work 
in this direction is promising and can be incorporated 
seamlessly in our framework. In addition, there are 
uncertainties in the equation of state, and in chemical 
and thermal freeze-out assumptions and parameters. 
We have not attempted a fine tuning of parameters - 
the sensitivity of our results to various parameters will 
be addressed in a follow up work. Despite these caveats, 
the successful description of a wide range of data in our 
model provides a framework to nail down key aspects of 
the complex dynamics of heavy ion collisions. 
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